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Compressive Sensing of Large-Scale Images: 
An Assumption-Free Approach 

Wei-Jie Liang, Gang-Xuan Lin, and Chun-Shien Lu 


Abstract —Cost-efficient compressive sensing of big media data 
with fast reconstructed high-quality results is very challenging. 
In this paper, we propose a new large-scale image compressive 
sensing method, composed of operator-based strategy in the 
context of fixed point continuation method and weighted LASSO 
with tree structure sparsity pattern. The main characteristic 
of our method is free from any assumptions and restrictions. 
The feasibility of our method is verified via simulations and 
comparisons with state-of-the-art algorithms. 

Index Terms —Compressed sensing. Convex optimization, 
Large-scale Images, Sparsity. 

I. Introduction 

C OMPRESSIVE sensing (CS) ffl, E), 0 of sparse 
signals in achieving simultaneous data acquisition and 
compression has been extensively studied in the literature. In 
the context of CS, we usually let u denote a /c-sparse signal of 
length n to be sensed, let $ of dimensionality mxn represent 
a sampling matrix, and let y be the measurement of length 
m, where k < m < n and 0 < — < 1 is defined as the 
measurement rate. At the encoder, random projection, defined 
as: 

y = <^u, (1) 

is conducted on the original signal u via $ to obtain the 
measurement vector y. At the decoder, u can be recovered 
based on its sparsity by means of convex optimization or 
greedy algorithms. 

Compressive sensing has been widely studied for ID signals 
and 2D images with reasonable sizes. However, efficient 
compressive sensing of big media data without making any 
assumptions has not been found in the literature. We can 
foresee that when signal length n is large, almost all existing 
CS algorithms will run out of memory and/or the computations 
will be overloaded. To address the need of compressive sensing 
of big media data in both efficient sensing and recovery 
aspects, we study a new scheme for big media signals. Please 
note that the storage consumption for a random matrix is 
unacceptable because it needs 2G memory to store a random 
matrix with (small) size n = (128 x 128)^. 

A. Related Work 

In a, Duarte and Baraniuk introduce Kronecker product 
to model multidimensional compressive sensing of signals 
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and propose a Kronecker compressive sensing (KCS) method. 
They prove that the Kronecker sensing matrix and Kronecker 
dictionary possess mutual incoherence property (MIP) and 
restricted isometric property (RIP). Nevertheless, the practical 
use of KCS is greatly prohibited because the vectorization of 
multidimensional signals and the use of joint sensing involve 
a (very) large Kronecker product-based sensing matrix. 

In 0, a multiway compressive sensing (MWCS) method 
for sparse and low-rank tensors is proposed. Although MWCS 
achieves more efficient reconstruction, its performance relies 
heavily on tensor rank estimation, which is NP-hard. A gener¬ 
alized tensor compressive sensing (GTCS) method for higher- 
order tensors has been proposed in IQ. GTCS is demonstrated 
to be comparable to KCS in recovery accuracy and be greatly 
faster than KCS in recovery speed. 

While the previous studies consider some tensor operations 
like Kronecker product, CP decomposition, and Tucker model 
within the framework of compressive sensing, the sparsity 
pattern inherent in the big media data/tensor (like 2D image 
and 3D video) has not been fully explored. Recently, Caiafa 
and Cichocki I?) exploit Kronecker product and block spar¬ 
sity to develop a so-called N-BOMP (N-way block OMP). 
However, we find, as also indicated in subsection 7.2.1 of Q, 
that for a 2D image it is pre-processed in advance to possess 
perfect block sparsity pattern in that the important/significant 
coefficients in some transform domain fall within the specified 
block sparsity pattern while other insignificant coefficients are 
entirely removed. Under the situation, N-BOMP is able to 
obtain reconstruction quality far better than the existing tensor 
CS algorithms. 

Recently, Caiafa and Cichocki IS present a fast tensor 
compressive sensing method, which no longer assumes certain 
sparsity pattern and does not involve iterations, making it 
suitable for large-scale problems. However, it assumes that 
the signal to be sensed and recovered has low multilinear-rank, 
leading to redundant sensing, which means that under the same 
measurement rate the reconstructed quality is (remarkably) 
lower than other CS algorithms. 

In a, we previously propose the use of tree-structure spar¬ 
sity pattern (TSSP) in tensor compressive sensing. TSSP can 
help to fast find significant wavelet coefficients. Its weakness 
is that there does not exist a fast recovery algorithm that can 
exploit TSSP. In this paper, we conquer this problem. 

Hale et al. m derive the optimality conditions of convex 
optimization problem as a fixed-point equation. Later, Wen et 
al. ifTTllfT^ provide a fast algorithm, called PPC_AS, to solve 
the convex optimization problem. Specifically, their sensing 
matrix for convex optimization problem is chosed as a partial 
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DCT matrix, whereas the fast algorithm is based on active-set 
strategy, which can solve the fixed-point equation in large- 
scale problems. 


B. Motivation and Our Contributions 

As indicated in Eq. when the signal length becomes 
large enough, the storage and computation of using the sensing 
matrix become an obstacle. To deal with sensing and recovery 
of large-sized images, we do not follow the convention of di¬ 
viding a large image into several small blocks ca, na, ca, 
m, wherein each small block can be separately processed. 
This will incur either blocky effects in recovery or required 
calibration in block-based sensing. 

Based on the above concerns, we explore the sensing 
strategy ifTTl and recovery strategy ifTOl . ifTTl . ifT^ . ifTSl that 
can be operated in operators to speed up computation and save 
storage without needing any assumptions/restrictions. In this 
paper, the recovery strategy is based on exploiting fixed-point 
equation Col to solve the convex optimization problem instead 
of greedy algorithms that rely on matrix computation, which 
is prohibited in the context of big data compressive sensing. 
Specifically, the recovery method is based on the FPC_AS 
algorithm CD, d that solves the fixed-point equation. To 
recover the large-scale data, the authors in d, CD, d 
propose an active-set algorithm and Milzarek et al. d 
propose a globalized semismooth Newton method to solve 
large-scale ^i-norm optimization problems, where partial DCT 
matrix is adopted as the sensing matrix for fast sensing. 
However, the signal to be sensed in d is assumed to be 
sparse in the time domain (that is the sparsifying basis is the 
identity matrix). 

In this paper, we propose to use a kind of random Gaussian- 
like matrix, called Structurally Random Matrix (SRM) d, as 
the sensing matrix, and wavelet as the sparsifying basis to deal 
with large-scale images. Our choice can satisfy RIP and MIP 
in CS. Basically, our method can be viewed as an extension 
of ifTSll to compressive sensing of large-scale signals that are 
not sparse themselves but such extension is not trivial at all. 

In addition to the above, for sparse recovery of big images, 
the sparsity pattern plays an important role. With an eye on 
the natural characteristic of tree-structure relationship among 
wavelet coefficients that are popularly used to represent media 
data, we propose to explore tree-structure sparsity pattern 
(TSSP) in big data compressive sensing. When TSSP is 
considered (with the way different from ||9l), our strategy is to 
assign smaller weights to the wavelet coefficients at the lower 
frequencies and larger weights to those at the higher frequen¬ 
cies under the framework of convex optimization. Specifically, 
we explore a weighted-LASSO algorithm for sparse recovery 
of big images from sensed measurements. 


II. BIG IMAGE COMPRESSIVE SENSING 

In this section, we describe the proposed method, wherein 
fixed point strategy is adopted to solve the Lasso problem and 
the step size is controlled by quasi-Armijo rule. 


A. System Model 

Let U G be the 2-dimensional image. It can be 

sparsely represented via certain dictionaries T'l and T '2 as: 

[/ = ( 2 ) 

where X is sparse with respect to T'l and 4*2. Here, we 

reshape 2D signal to ID vector. Based on Eq. Q, we have 

y = <i)M = ‘hT'x = Ax, (3) 

where A = G n = N^, sparsifying basis T* = 

'^2 C) 'ki G u = vec{U), and x = vec{X) G M". 

The main problem is to reconstruct the vector x from the 
measurement y in Eq. (0. A well-established approach for 
the reconstruction is an optimization method which we call 
the LASSO (or penalized least-squares) problem: 

i||y-Ax||2, (4) 

Z 

where A > 0 specifies the penalty of sparse level of x. To ease 
discussion later, we set J^(x) = A ||x||^ + ^Wu — ^ 2 :|| 2 . 

B. Sensing Matrix Design 

The common use of random Gaussian matrix as the sensing 
matrix leads to problems of storage and computation cost. 
Although storage consumption can be overcome by using a 
seed to generate random Gaussian, it still encounters high 
computational cost. Nguyen et al. in CD propose a frame¬ 
work, called Structurally Random Matrix (SRM), with: 

$ = DFR, (5) 

where D G is a sampling matrix, F G is an 

orthonormal matrix, and R G is a uniform random 

permutation matrix. Since the distributions between a random 
Gaussian matrix and SRM’s $ are verified to be similar, we 
choose Eq. 0 as the sensing matrix for our use. 

In this paper, we set F to the Discrete Cosine Transform 
(DCT) due to its fast computation and cost-effectiveness. 


C. Fixed Point Method with Quasi-Armijo Rule 

Due to the convexity of the function F(x) in Eq. Q, the 
global minimum solution is exactly the critical point of F{x). 
In m, the authors derive that the critical point of F{x) 
exactly belongs to the solution set of fixed point equation as 


J = 


X = Sr\{Gr{x)) 


( 6 ) 


where r > 0 is arbitrarily fixed, 

r Gr(x) = X — tA'^{Ax — y), 

\ ^ tAjtA](^), 


and = min {max {—c, f} , c} is the projection onto 

the interval [—c, c]. Then the fixed point iteration is defined 
as: 


xfe+i _ SrxiGrix^)) with T > 0. (8) 


Let x^ be the current iterate and let — x^ be a 

direction that is generated by Eq. Then we can calculate 
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2 ,fc+i i^y where ak is the step size controlled by 

a quasi-Armijo rule. The algorithm is depicted in Algorithm 
[2 where we extend it to deal with big images that they 
themselves are not sparse. 


Algorithm 1 ifTSl Fixed point method with quasi-Armijo rule. 
Input: The initial iterative point, = 0 S K”; The 

shrinkage parameter, r > 0; The quasi Armijo’s step size 
parameter, (3 € (0,1); A constant, 7 G (0,1); The 
weighted parameter of LASSO, A > 0; An initial 
iterative step, k = O', 

Output: The k'^ iterative point, 

1 : Calculate the direction: = Srx{Gr{x’^)) — x^', 

2: while ^ 0 do 
3 : Calculate = 

{A^ {Ax - 2/)) -f A (Gr (a;''))||i - ||a:''||J; 

4 : Choose a maximal quasi-Armijo step size 

ak & { 1 , /3, /3^, ■ •.} such that 
Tix’^ + akd'^) - Tix’^) < akjAsk; 

5 : = x'^ + akd!^', 

6: end while 
7 : return x^; 


D. Tree Structure Sparsity in Convex Optimization 


For improving the quality of reconstruction, we refer to an 
iterative reweighted Zi-norm minimization (IRWLl), which is 
proposed in ESI. For a given diagonal weighted matrix W G 
Ijnxra, jjjg convex optimization problem in Eq. Q can be 
relaxed with x = Wx and A = as: 

2 


min A I 




Ax — y 


(9) 


The main difference between our model in Eq. (|^ and IRWLl 
is that the weighted matrix W is determined by tree structure 
sparsity pattern (TSSP). 

More specifically, TSSP is yielded by separating the wavelet 
coefficient for 2D image into support and not-support sets. We 
adopt the wavelet, as provided in the source code of Q, as 
the dictionary T* and apply it to an image with S levels to 
obtain the subbands {LLs, LHs, HLs, HHs, LHs-i, ■ ■ ■, HHA, 
where L and H denote low and high frequencies, respectively. 
The resultant wavelet coefficients are weighted according to 
which levels of subbands they are located in. 

Algorithm [^describes how to solve Eq. (j^. Eirst, since the 
coefficients in LLs are significant, the corresponding indices 
in W are all reserved and set to 0.1 (as in initialization part 
of Algorithm S, and other diagonal elements are set to 1. 
We solve Eq. Q (Step in Algorithm to yield the initial 
solution. 

Second, we check the entries from the subbands LHs, HLs 
and HHs that could be the roots of evolving trees (Steps 
andj^in Algorithm]^. They will be put in the queue Q if they 
are supports (corresponding coefficients are large enough (Step 
1^ in Algorithm |^). Eor each element in Q, if it is checked 
to be a support, its children will be put in Q (Steps and 
in Algorithm]^. This process is repeated until Q is empty to 
complete the generation of TSSR 


To construct W, its diagonal entries are decided by TSSR 
We expect that the wavelet coefficients, solved by Eq. 
are located on the tree decided by TSSR. Since the decision 
variable with small weight will derive the large wavelet 
coefficient, together with the fact that the energy is decreasing 
from level S to level 1, the weights are empirically set to: 

dmg(W)^ = 0A(S-s + l) if i € Is, (10) 

where A is the subset of LHg U HLgU HHs and is decided 
by TSSR. 

Einally, we solve Eq. (0 with weighted matrix W in Eq. 
([T0|) to obtain the final solution (Stepin Algorithm]^. 


Algorithm 2 IRWLl with TSSR. 


Input: The initial weighted iterative point, a; = 0 G M"; The 
initial truncated iterative point, i = 0 G M”; The 
weighted matrix, W = Inxn', The index set of LLs, 
Is+i', The empty sets. Is, f ^ s < S', The percentage, 
p%; A tolerance, e > 0; 

Output: The final output, x 
1 : 

2 : Solve Eq. by Alg. [T] with weighted matrix W to 
obtain the solution a:*; 


Set the weighted of LLs as diag{W)\j^^^ = 0.1; 


7 : 

8 : 

9 : 

10 

11 

12 

13 

14 

15 
16 : 

17 : 


Set the truncated variable a; L = a;* L ; 

Jjs+l IJS+l 

Calculate residue r = y — Ax; 

Calculate correlation = A^ r , i = 1,2,... ,n; 

Let Is collects the indices with the first p% largest 
correlations in LHs U HLs U HHs. Set 
diag{W)\^ =0.2; 
for s = S' — 1, S' — 2,..., 1 do 

Construct A as the children of Is+i, set diag{W)\^ 
according to Eq. 

Solve Eq. (j^ by Alg. 1 to obtain the solution x*; 
for j G Is do 

if I a;* I < e then 

remove index j and set diag{W)j = 1; 

end if 
end for 
end for 

Solve Eq. by Alg. 1 with final weighted matrix W 
and output the solution x*; 

return; 


E. Memory Cost and Computational Complexity 

The main computational cost of solving Eq. (j^ comes from 
the computation of matrix-vector multiplications: 


Ax = <l>'lfW-^x = DFRWW- 


including Gr{x) = x — {^Ax — y^ and J" (a;*) = 


A| 


y — Ax 


, where D, F, and R are defined in Sec. 


|II-B[ W is an inverse wavelet transform, and W is defined in 
Sec. IlLDl 

Since the size of a: is n = the storage cost for the 
matrices F, R, W, and W is n? = N'^, and is to x n for matrix 
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D. For example, if the image we aim to reconstruct is of size 
128 X 128, then the matrices, F, R, W, W, cost around 8GB 
memory in total. However, if we consider that the permutation 
matrix R and diagonal matrix W can be represented by n = 
entries, around 4GB are enough. As a result, the limited 
storage leads to the restriction of image size. 

In order to speed up the computation of linear transforma¬ 
tion (i.e., matrix-vector multiplications here), we resort to lin¬ 
ear operator in MATLAB or reformulation, as described below. 
In other words, the memories required to store the matrices, 
mentioned above, can be remarkably reduced accordingly such 
that our method can be adaptive to large images. 

• Weighted matrix W G 

Wx = w o X, where w = vec{W) and o denotes the 
Hadamard product. The memory cost of W is defined as 
costM{W) = n due to Hadamard operation. 

• Inverse wavelet matrix W G 

Wx = (W2 ® Wi) a; = t;ec(WiXW|’). 

We can see that costM{W) = 2n due to the use of 
Kronecker product. 

• Random permutation matrix R G 

Rx = TZ{x), where the operator TZ{-) randomly permutes 
the indices of vector x. Thus, costM{R) = n. 

• Discrete Cosine Transform (DCT) F G 

F{x) = dct{x), where DCT can be calculated by the 
dct operator in MATLAB, which is speeded up by Fast 
Fourier Transform. Thus, we have costM{F) = n. 

• Partial random permutation D G 

Dx = ^{x), where the operator !?(•) randomly 

chooses m indices from n components in vector x. So 
costMiF) = m. 

Therefore, the memory cost of our method is in total m + 6n. 

On the other hand, the computational complexity for cal¬ 
culating Ax by matrix-vector multiplications and by linear 
operator are compared as follows. 

• Matrix-vector multiplication: Since 

Ax = DFRWW~^x = ^yVW~^x, the time complexity 
for individual matrix multiplication is: 


matrix 

$ = DFR 

W 

w-^ 

complexity 

0(mn) 

0{n^) 

0{nA 


• Linear operations or reformulations: The time complexity 


for individual operation is: 


operation 

D 

F 

R 

W 

W-^ 

complexity 

0{m) 

0 (nlog(n)) 

0{n) 

O(n^F) 

0(n) 


We can see that the computational complexity 0{n^) of 
matrix-vector multiplications is reduced to 0{rA^^) of linear 
operations. 


E Convergence Analysis 


In this section, we study the convergence of the solution 
sequence {x„}, which is generated by AlgorithmBased on 
llQl, we choose 

T G (0, 2/A„iax) , 

which guarantees that both two functions in Eq. 0 are 
nonexpansive, where 


Ama,v .— max An 


\y-Ax\\ 


Theorem 1. Let {xn\ be a sequence generated by Algorithm 
Assume that J Then {xn} converges to a point in J. 

The strategy of proving Algorithm is to prove that the 
sequence {xn} is Cauchy, with the fact that every Cauchy 
sequence converges in complete space M". 

By Theorem the sequence {x„} generated by Algorithm 
converges to a fixed point x of the fixed point equation Eq. 

i.e., the optimal solution of Eq. 0 . 

Moreover, Algorithm is designed based on Algorithm [T] 
with the weighted matrix constructed in terms of tree structure, 
that is. Algorithm is conducted by iteratively performing 
Algorithm [T] S'— 1 times. Thus, the convergence of Algorithm 
is guaranteed by the convergence of Algorithm We 

show the normalized function errors ^ vs. number 

__ '' (^21 ) 

of iterations in Eig. [H where the measurement rates range 
from 10% to 30%. We can see that the algorithm converges 
more and more fast (with less number of iterations) as the 
measurement rates increase. 

III. Simulation Results 

All simulations were conducted on PC equipped with Win¬ 
dows 7 with 3.4 GHz Intel Core i7 CPU and 8 GB RAM. We 
compare the proposed method, N-BOMP j7], and 0 in terms 
of CPU time for CS recovery and reconstruction quality. Here, 
we take 2D image of large sizes (ranging from 1024 x 1024 to 
4096 X 4096) for compressive sensing and recovery, and show 
some results for images with different degree of sparsity in this 
section. The images are shown in Pig. [T] (a). They are (from 
left to right) Paint, Man (the first two adopted in fT)), airport 
(a commonly used image in the literature), and a synthetic 
aperture radar (SAR) image (a TerraSAR-X image downloaded 
from http://www.geo-airbusds.com/en/23-sample-imagery). 

A. Parameter Setting 

The sensing matrix was a Structurally Random Matrix HU, 
the Daubechies wavelet was used as the sparsifying basis, and 
fixed point method with quasi-Armijo rule was adopted for 
CS recovery. Note that the images to be sensed were not pre- 
processed in advance. 

B. Performance Comparison 

As mentioned previously, the four images of size 1024 x 
1024 shown in Pig. [T](a) were used as the targets for sensing. 
Several measurement rates (MRs), ranging from 5% to 50%, 
were selected for simulations. The CPU time, the reconstruc¬ 
tion quality measured in Peak-Signal-to-Noise-Ratio (PSNR), 
and the reconstruction quality measured in structural similarity 
(SSIM) EOll were adopted as the criteria for comparison 
among our method (Algorithm^, N-BOMP iTl, and IH. Some 
results are shown in the following figures. 

Prom Pig. m to Pig. 0 we can find that under the same 
measurement rates, our method obtains the best reconstruction 
qualities in term of PSNR and SSIM within reasonable recov¬ 
ery time. Although our method needs more computational time 
than the other two methods, it does not rely on any impractical 
assumptions and restrictions (block sparsity or low rank) and 
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(al) (a2) (a3) (a4) 

iteration vs. relative error iteration vs. relative error iteration vs. relative error iteration vs. relative error 



(bl) (b2) (b3) (b4) 

iteration vs. relative error iteration vs. relative error iteration vs. relative error iteration vs. relative error 



(cl) (c2) (c3) (c4) 

iteration vs. relative error iteration vs. relative error iteration vs. relative error iteration vs. relative error 




(dl) (d2) (d3) (d4) 


Fig. 1. The comparison of normalized function error — ^ —- (r/-axis) versus number of iteration (a:-axis) under different measurement rates in row 

J- (a; ) 

(b) 10%; row (c) 20%; and row (d) 30% for the corresponding four images in row (a). 


is flexible in practice. In particular, we do not think it is 
meaningful and impressive to quickly yield inaccurate results. 

Furthermore, for visual comparison, we actually observe 
that the image details can be well recovered by our method. 
This also explains the usefulness of tree-structure sparsity 
pattern imposed in the ti-norm minimization formulation (Eq. 
Q). In fact, our simulations also show that when weighting 
is imposed, the PSNR gain of 0.5 ~ 1 dB can be obtained, 
when compared to its non-weighting counterpart. 

Finally, when the popularly used CVX software is consid¬ 
ered for CS recovery, our results show that basically CVX 
cannot deal with images larger than 100 x 100 due to high 
cost of memory and computation. 


IV. Conclusion 

Compressive sensing of large-scale images is remarkably 
challenging due to the constraints of storage and compu¬ 
tation. In this paper, we propose a new method of large- 
scale image compressive sensing based on exploring a fixed- 
point weighted-LASSO algorithm without depending on any 
assumption or preprocessing of sparsity pattern in images. 
Convergence analysis is also provided to confirm the conver¬ 
gence of our iterative scheme. 

Appendix A 
Proof of Theorem 1 

To prove Theorem 1, we need the following Lemmata. 




















































6 



Fig. 2. Comparison of time/PSNR/SSIM of the “Paint” image among 
Algorithmic NBOMP Q, and H). 


Fig. 4. Comparison of time/PSNR/SSIM of the “Airport” image among 
Algorithmic NBOMP Q, and g). 



Fig. 3. Comparison of time/PSNR/SSIM of the “Man” image among 
Algorithmic NBOMP (T), and g). 

Lemma 1. The operator in Eq. ^ is nonexpansive 

(Lemma 3.2 in mi); Gr{-) is nonexpansive for an appropriate 
parameter r (Lemma 4.1 and Eq. (4.3) in m); and thus 
o Gr) (•) is as well. Moreover, o Gr) (•) is contin¬ 
uous. 

Lemma 2. In a normed vector space {X, |j • |j), the following 
identical equation is called the parallelogram law: 

2\\xr + 2\\yr = \\x + yr + \\x-y\\^ Vx, y S X 

Lemma 3. Let {a^},{5^} be non-negative sequences, and 
{a^} does not converge to 0. If lim = 0, then 

k—^oo 

lim inf = 0. 

k—foo 


Fig. 5. Comparison of time/PSNR/SSIM of an SAR image among Algorithm 
1C NBOMP 0), and g). 


Proof: Since {a^} does not converge to 0, there exists 
ei > 0 s.t. for all j S N, we can find Uj > j so that |a"j —0| > 
ei- 

We have 

0 = lim > lim > 0, 

j—foo j—foo 

which implies lim b^^ = 0. In other words, the limit of 

j-t-oo 

subsequence {&”j } exists, and hence 

lim inf 6"^ = lim 6"^. 

k—foonj'>k j—>-oo 
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By the fact that {b^}j>k 3 {b'^^}n^>k, we have 


inf b> < inf 6 ”’. 

j^k '^j'^k 


Thus 


0 < liminf = lim inf V < lim inf 6 ”^ = lim 6 "^ = 0. 

k—>-oo k—^ooj^k k—^oorij^k j—>-oo 

Finally, we have lim inf = 0. ■ 

A;—>-oo 

In the following, we prove Theorem 

Proof: First we show that the limit lim \\x^ — z\\ exists. 

k—foo 

Let z G and define = Sr^i o Gt{x^)- Then 
x^+^ - z\\^ 

x^ + ak {V^ - x’^) - z\ 

(1 - (Jk)x'^ + (TkV'^ - z [2 
(1 - (Tk){x^ - z) + (Tfe {'P'" - z)\\^ 

< (1 - CTfc) IlcCfc - ZII 2 + CTfc \\V^ - zIL 

< (I - a^)\\x’^ - z\\^ + ak\\x^ - z\\^ 

(since z G J and by Lemma 


= \\x — Z I 


llcc'=-z|| 


is monotone decreasing and 


Hence, the sequence 
lim \\x^ — z\\ exists. 

k—^oo 

Second we show that the sequence is Cauchy. By 

Lemma 1^ we have 


||a:^-x '=||2 

= ll(x^-z)- (x'=-z )||2 

= 2||x^ - z||2 + 2||x*^' - zf - \\x^ + x'= - 2 z||2. 
Let lim ||x^ — z|| = c. We have 

k—¥oo 

\\\x^+x’^ - 2 z|| - 2 c| 

= |||(:,^-^) + (x'=-z)||-c-c| 

< |||x^-d|-d + |||x'=-z||-d 


( 11 ) 


and 


lim 

h,k—f(x> 


lx'* - 2 z|| - 2 c| 


< lim nix'* — z|| — c| + lim |||x'* — z|| — c| 

h—foo k—foo 

= 0, 


which means 


lim ||x'* + x'* — 2 z|| = 2 c. 

h.k—foo 


( 12 ) 


By Eq. ([nj and Eq. ( [T^ , we obtain 

lim II x'* — x ''||2 = 2 c^ + 2 c^ — ( 2 c)^ = 0 , 

h,k—foo 

and hence the sequence {x'*} is Cauchy. Therefore, {x'^} is a 
convergent sequence. 

Third we prove that limit x* = lim x^ G J. Since x^^^ = 

k—^oo 

x^ + <7k{P^ — x^), we have 

(Tfc||7^"--x'=|| = ||x"-+i-x'=||, 

and then 

lim CTfe IIT^'^ — x'* II = lim ||x''“'''-— x'^ || 


(13) 


k—foo 


k—¥oo ' 


lim — x^) 

k—foo 


X — X\\ 


= 0. 


By the fact that sequence {ok} is decided by Stepfor each 
iteration in Algorithm[T] each ak is mutually independent, and, 
thus, {ak} does not converge. By Lemma we have 

liminf||lP'*-x'*|| =0. (14) 

k—^oo 


Next we aim to show that lim |j7^'* — x'^H exists. Since 

k^oo 


< 


-pk+l _ j.fc+1 II 

1P'*+1 - x'* - crfc(lP'= - x'=)|| 

- {I - ak)x^ - akV^W 

'pk+l _'pk 

jyk+l _ pfc (1 _ cTfeX'P'* - x'=)[| 
pk+1 _ j,k\\ ii^fc - a;H| , 


by the fact that S-r^i o Gr is nonexpansive, that is 
||pfc+l _'pk\\ < ll^-fc+l _2.fc|| ^ 


According to Eq. ( [T3] l, it follows that 

||pfc+i_pfc||_^(l_^^)||pfc_^fe 

< ||x'=+i -x'*'|| +(l-tTfc)||lP''-a:'=|| 
= (Tfc||lP'=-x'=|| +(l-afc)||lP'*-x'*|| 
= || 1 P'= -x'*||. 


This implies that {||7^'' — x'*||} is a decreasing sequence and 
bounded below by 0. Therefore, the sequence {||7^'^ — x'^||} is 
convergent and Eq. leads to 

lim lllP'^-x'^ll = 0 . 

k—^oo 

Einally, according to Lemma both functions S-r^oGr and 
II'll are continuous. We have 


0 = lim 117^'= -x'* I 

k—foo 


lim (V^ - x'=) 

k—foo 


( lim x'^ j 

— lim X* 

\k—foo J 

k—^oo 


Sr^Gr{x)-X*\\. 


Hence x* = lim x^ is a fixed point of function S'^-^ o Gr, 

k—^oo 

and we complete the proof. ■ 
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